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1.  Abstract 

This  project  is  about  the  design,  analysis  and  application  of  high  order  accurate 
and  nonlinearly  stable  finite  difference  (including  finite  volume),  finite  element  and 
spectral  algorithms  for  computing  solutions  of  partial  differential  equations  which 
are  either  discontinuous  or  with  sharp  gradients.  Algorithm  development,  theoretical 
study  about  stability  and  convergence  of  the  algorithms,  investigation  about  efficient 
implementation  including  parallel  implementations,  and  applications  in  compressible 
and  incompressible  gas  dynamics  and  in  semiconductor  device  simulations,  are  per¬ 
formed.  The  achievement  strengthens  our  objective  to  obtain  powerful  and  reliable 
high  order  numerical  algorithms  and  use  them  to  solve  problems  containing  discon¬ 
tinuous  solutions,  especially  those  of  army  interest. 


2.  Statement  of  the  Problem  Studied 

The  problems  studied  in  this  project  involve  numerical  solutions  of  partial  dif¬ 
ferential  equations,  mainly  hyperbolic  type  or  convection  dominated  parabolic  type 
equations,  with  solutions  which  are  either  discontinuous,  or  with  discontinuous  deriva¬ 
tives,  or  containing  sharp  gradient  regions  which  are  difficult  to  be  completely  resolved 
on  today’s  computer.  The  methods  we  investigate  fall  into  the  category  of  “shock 
capturing”  schemes,  which  means  that  these  methods  try  to  capture  shocks  or  other 
types  of  discontinuities  and/or  sharp  gradient  regions  with  a  relatively  coarse  grid, 
rather  than  completely  resolving  them.  These  methods  are  useful  when  it  is  either 
impossible  or  too  costly  to  completely  resolve  certain  solution  structure.  High  order 
accurate  finite  difference,  finite  element  and  spectral  methods  have  all  been  consid¬ 
ered. 

Our  approach  is  to  explore  the  robustness  and  efficiency  of  high  order  numerical  al¬ 
gorithms  for  nonsmooth  problems  both  through  theoretical  guidance,  often  obtained 
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with  rigorous  proofs  on  simplified  model  problems,  and  through  numerical  experi¬ 
ments  on  real  application  problems.  We  do  not  try  to  modify  algorithms  just  for 
the  purpose  of  convergence  proofs,  if  such  modifications  are  not  justified  by  numer¬ 
ical  experiments.  For  finite  difference  schemes,  we  are  exploring  the  very  efficient 
ENO  schemes  of  Shu  and  Osher  based  on  point  values,  numerical  fluxes,  and  nonlin- 
early  stable  high  order  Runge-Kutta  time  discretizations,  and  the  recently  developed 
WENO  (Weighted  ENO)  schemes.  For  finite  element  methods,  we  are  exploring  the 
Runge-Kutta  discontinuous  Galerkin  methods  of  Cockburn  and  Shu,  which  combine 
the  advantage  of  finite  elements  (weak  formulation,  automatic  energy  stability,  easy 
handling  of  complicated  geometry  and  boundary  conditions)  with  features  of  high 
resolution  finite  difference  schemes  (approximate  Riemann  solvers,  limiters).  Effec¬ 
tive  ways  to  handle  viscous  terms  are  being  investigated.  For  spectral  methods,  we 
are  exploring  reconstruction  techniques  of  Gottlieb  and  Shu  to  apply  spectral  approx¬ 
imations  to  discontinuous  functions  and  still  obtain  uniform  spectral  accuracy. 

We  have  been  cooperating  with  Dr.  Rupak  Biswas  of  RIACS  and  Dr.  Roger 
Strawn  of  US  Army  AFDD,  at  NASA  Ames  Research  Center,  on  the  investigation 
of  developing  high  order  high  resolution  numerical  methods  for  the  simulation  of 
helicopter  rotor  blade  motion.  This  is  a  very  demanding  simulation  since  the  vortex 
motion  after  interaction  with  shocks  should  be  sustained  for  long  time/distance,  and 
only  high  order  methods  with  low  dissipation  can  achieve  this.  We  have  already 
performed  simulations  for  a  model  problem,  namely  a  moving  vortex  in  a  uniform 
flow  which  is  an  exact  solution  of  the  Euler  equations.  Both  2D  and  3D  simulations 
are  performed,  with  vortex  moving  in  different  directions  and  with  both  uniform  and 
nonuniform  (but  smooth)  Cartesian  grids.  The  results  clearly  show  the  advantage  of 
using  high  order  methods  for  this  model  problem.  We  are  currently  in  the  process 
of  combining  the  WENO  scheme  (as  the  background  solver)  and  the  existing  inner 
solver  near  the  blade  which  is  of  finite  volume  type. 


3.  Summary  of  Research  Results 

Research  has  been  performed  in  all  areas  listed  in  the  original  proposal,  and 
progress  and  results  consistent  with  the  original  objectives  have  been  obtained.  There 
are  28  publications  (among  them  22  in  refereed  journals)  resulting  from  this  project, 
see  Section  4  for  a  list  of  them. 

About  high  order  essentially  non-oscillatory  (ENO)  finite  difference  schemes,  jointly 
with  Yanni  Zeng,  we  have  applied  ENO  method  to  the  viscoelastic  model  with  fad¬ 
ing  memory  [15]  (all  the  numbering  of  references  are  according  to  that  of  Section  4). 
The  memory  term  is  treated  by  introducing  new  variables  and  rewrite  the  system  by 
adding  more  differential  equations  but  without  explicit  memory  terms.  The  appear¬ 
ance  of  the  memory  terms  regularizes  the  solution  somewhat,  and  in  many  cases  it 
is  still  a  theoretically  open  question  whether  shocks  will  develop  from  smooth  initial 
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data.  We  have  performed  theoretical  analysis  about  the  linearized  system  for  large 
time,  and  have  applied  ENO  scheme  to  study  the  nonlinear  system  for  both  local  time 
and  large  time.  The  high  order  accuracy  and  sharp,  non-oscillatory  shock  transition 
allow  us  to  obtain  fine  resolution  for  tens  of  thousands  of  time  steps,  and  to  study 
the  shock  interactions  after  the  formation  of  shocks. 

In  [4],  which  results  from  an  invited  talk  in  the  Third  International  Colloquium 
of  Numerical  Analysis,  we  have  compared  ENO  finite  difference,  finite  volume  and 
discontinuous  Galerkin  finite  element  methods,  in  terms  of  their  shock  resolution,  grid 
orientation  effects,  cost,  and  ability  to  handle  complicated  geometry  and  boundary 
conditions. 

Also  on  finite  difference,  Harabetian,  Osher  and  Shu  have  investigated  a  novel 
Eulerian  approach  for  simulating  vortex  motion  using  a  level  set  regularization  proce¬ 
dure  [12].  Our  approach  uses  a  decomposition  of  the  vorticity  of  the  form  £  =  P{<p)r], 
in  which  both  </?  (the  level  set  function)  and  rj  (the  vorticity  strength  vector)  are 
smooth.  We  derive  coupled  equations  for  (p  and  77  which  give  a  regularization  of  the 
problem.  The  regularization  is  topological  and  is  automatically  accomplished  through 
the  use  of  numerical  schemes  whose  viscosity  shrinks  to  zero  with  grid  size.  There  is 
no  need  for  explicit  filtering,  even  when  singularities  appear  in  the  front.  The  method 
also  has  the  advantage  of  automatically  allowing  topological  changes  such  as  merging 
of  surfaces.  Numerical  examples  including  two  and  three  dimensional  vortex  sheets, 
two  dimensional  vortex  dipole  sheets  and  point  vortices,  are  given.  To  our  knowledge, 
this  is  the  first  three  dimensional  vortex  sheet  calculation  in  which  the  sheet  evolution 
feeds  back  to  the  calculation  of  the  fluid  velocity. 

Application  of  ENO  scheme  to  the  study  of  shock  longitudinal  vortex  interaction 
problem  is  carried  out,  jointly  with  Erlebacher  and  Hussaini  [16].  We  have  studied 
the  shock/longitudinal  vortex  interaction  problem  in  axisymmetric  geometry.  Linear 
analysis,  shock  fitting  code,  and  shock  capturing  ENO  are  used  in  different  parameter 
range,  to  study  various  cases  of  nearly  linear  regime,  weakly  nonlinear  regime,  and 
strong  nonlinear  regime.  Vortex  breakdown  as  a  function  of  Mach  number  ranging 
from  1.3  to  10  is  studied,  extending  the  range  of  existing  results.  For  vortex  strengths 
above  a  critical  value,  a  triple  point  forms  on  the  shock,  leading  to  a  Mach  disk.  This 
leads  to  a  strong  recirculating  region  downstream  of  the  shock.  It  is  found  out  that  a 
secondary  shock  forms,  to  provide  the  necessary  deceleration  so  that  the  fluid  velocity 
can  adjust  to  downstream  conditions  at  the  shock. 

The  important  issue  of  maintaining  positivity  of  density  and  pressure  for  high 
order  calculations  of  Euler  equations  of  compressible  gas  dynamics  is  considered  in 
[11]. 

As  an  important  development  about  the  high  order  essentially  non-oscillatory 
(ENO)  finite  difference  schemes,  Jiang  and  Shu  have  been  investigating  WENO 
(weighted  ENO)  schemes  [13].  WENO  is  a  modification  and  improvement  of  ENO 
schemes.  Instead  of  using  only  one  of  the  many  candidate  stencils  based  on  local 
smoothness  as  in  ENO,  WENO  uses  a  linear  combination  of  the  contribution  from 
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all  candidate  stencils,  each  with  suitable  nonlinear  weight.  The  weights  are  chosen 
so  that  in  smooth  regions,  they  are  close  to  an  optimal  linear  weight  which  gives  the 
highest  possible  order  of  accuracy  of  an  upwind-biased  linearly  stable  scheme.  Near 
shocks,  however,  those  stencils  which  contain  the  shock  are  assigned  weights  which  are 
essentially  zero.  Thus  WENO  resembles  a  linear  high  order  upwind  biased  scheme 
in  smooth  regions,  and  ENO  near  shocks,  with  a  smooth  numerical  flux  function. 
Another  advantage  of  WENO,  due  to  its  smoothness  of  fluxes,  is  that  convergence 
for  smooth  solutions  can  be  proven.  Also,  convergence  towards  steady  states  is  eas¬ 
ier  than  ENO.  Numerical  experiments,  including  shock  vortex  interaction  problems 
relevant  to  the  joint  work  of  the  PI  and  ARO  scientists  Drs.  Rupak  and  Strawn  on 
the  simulation  of  helicopter  rotor  blade  motion,  have  been  performed  with  very  good 
results. 

In  [26],  jointly  with  Siddiqi  and  Kimia,  we  have  developed  a  geometric  shock 
capturing  ENO  methodology  suitable  for  subpixel  interpolation  (subcell  resolution). 
The  method  gives  excellent  result  for  curve  evolution  and  representation  on  a  fixed, 
coarse  grid. 

In  the  lecture  note  [27],  we  describe  the  construction,  analysis,  and  applica¬ 
tion  of  ENO  (Essentially  Non-Oscillatory)  and  WENO  (Weighted  Essentially  Non- 
Oscillatory)  schemes  for  hyperbolic  conservation  laws  and  related  Hamilton- Jacobi 
equations.  ENO  and  WENO  schemes  are  high  order  accurate  finite  difference  schemes 
designed  for  problems  with  piecewise  smooth  solutions  containing  discontinuities.  The 
key  idea  lies  at  the  approximation  level,  where  a  nonlinear  adaptive  procedure  is  used 
to  automatically  choose  the  locally  smoothest  stencil,  hence  avoiding  crossing  dis¬ 
continuities  in  the  interpolation  procedure  as  much  as  possible.  ENO  and  WENO 
schemes  have  been  quite  successful  in  applications,  especially  for  problems  containing 
both  shocks  and  complicated  smooth  solution  structures,  such  as  compressible  turbu¬ 
lence  simulations  and  aeroacoustics.  This  lecture  note  is  basically  self-contained.  It 
is  our  hope  that  with  this  note  and  with  the  help  of  the  quoted  references,  the  reader 
can  understand  the  algorithms  and  code  them  up  for  applications.  Sample  codes  are 
also  available  from  the  author. 

About  the  discontinuous  Galerkin  methods,  Atkins  and  Shu  have  developed  a  dis¬ 
continuous  Galerkin  formulation  that  avoids  the  use  of  discrete  quadrature  formulas 
[21].  The  application  is  carried  out  for  one  and  two  dimensional  linear  and  nonlinear 
test  problems.  This  approach  requires  less  computational  time  and  storage  than  con¬ 
ventional  implementations  but  preserves  the  compactness  and  robustness  inherent  in 
the  discontinuous  Galerkin  method. 

In  [22],  jointly  with  Cockburn,  we  study  the  Local  Discontinuous  Galerkin  meth¬ 
ods  for  nonlinear,  time-dependent  convection-diffusion  systems.  These  methods  are 
an  extension  of  the  Runge-Kutta  Discontinuous  Galerkin  methods  for  purely  hyper¬ 
bolic  systems  to  convection-diffusion  systems  and  share  with  those  methods  their  high 
parallelizability,  their  high-order  formal  accuracy,  and  their  easy  handling  of  compli¬ 
cated  geometries,  for  convection  dominated  problems.  It  is  proven  that  for  scalar 
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equations,  the  Local  Discontinuous  Galerkin  methods  are  L2-stable  in  the  nonlinear 
case.  Moreover,  in  the  linear  case,  it  is  shown  that  if  polynomials  of  degree  k  are 
used,  the  methods  are  k- th  order  accurate  for  general  triangulations;  although  this 
order  of  convergence  is  suboptimal,  it  is  sharp  for  the  LDG  methods.  Preliminary 
numerical  examples  displaying  the  performance  of  the  method  are  shown. 

In  [25],  again  jointly  with  Cockburn,  we  extend  the  Runge-Kutta  discontinuous 
Galerkin  method  to  multidimensional  nonlinear  systems  of  conservation  laws.  The 
algorithms  are  described  and  discussed,  including  algorithm  formulation  and  practi¬ 
cal  implementation  issues  such  as  the  numerical  fluxes,  quadrature  rules,  degrees  of 
freedom,  and  the  slope  limiters,  both  in  the  triangular  and  the  rectangular  element 
cases.  Numerical  experiments  for  two  dimensional  Euler  equations  of  compressible 
gas  dynamics  are  presented  that  show  the  effect  of  the  (formal)  order  of  accuracy  and 
the  use  of  triangles  or  rectangles,  on  the  quality  of  the  approximation. 

Jointly  with  Hu,  we  present  a  discontinuous  Galerkin  finite  element  method  for 
solving  the  nonlinear  Hamilton- Jacobi  equations  [28].  This  method  is  based  on  the 
Runge-Kutta  discontinuous  Galerkin  finite  element  method  for  solving  coiiservation 
laws.  The  method  has  the  flexibility  of  treating  complicated  geometry  by  using  arbi¬ 
trary  triangulation,  can  achieve  high  order  accuracy  with  a  local,  compact  stencil,  and 
are  suited  for  efficient  parallel  implementation.  One  and  two  dimensional  numerical 
examples  are  given  to  illustrate  the  capability  of  the  method. 

For  spectral  methods,  jointly  with  David  Gottlieb,  we  have  been  continuing  on  our 
investigation  of  overcoming  Gibbs  phenomenon,  i.e.,  to  recover  exponential  accuracy 
from  a  spectral  partial  sum  of  a  piecewise  analytic  function.  We  have  investigated 
the  problem  of  recovering  exponential  accuracy  in  a  sub-interval,  with  the  assumption 
that  the  function  is  analytic  in  this  sub-interval  but  may  have  discontinuities  at  one 
or  both  boundary.  The  Fourier  and  Legendre  Galerkin  cases  are  considered  in  [10]. 
The  general  Gegenbauer  Galerkin  case  is  considered  in  [1].  The  collocation  case 
is  considered  in  [2].  We  review  the  history  of  the  Gibbs  phenomenon  as  well  as 
efforts  for  its  resolution  in  [17].  Also,  jointly  with  Peter  Wong,  we  have  performed  a 
detailed  numerical  study  about  numerical  accuracy  when  spectral  method  is  applied 
to  a  nonlinear  conservation  law  with  discontinuous  solutions  [5],  [14].  We  assess  the 
accuracy  of  Fourier  Galerkin  and  collocation  method  applied  to  Burgers  equation 
with  smooth  initial  condition  but  with  a  shock  developed  in  finite  time.  We  find 
that,  unlike  in  the  linear  PDE  case,  the  moments  with  respect  to  analytic  functions, 
in  particular  the  first  few  Fourier  coefficients,  are  no  longer  very  accurate.  However 
the  numerical  solution  does  contain  accurate  information  which  can  be  extracted  by 
a  post-processing  based  on  Gegenbauer  polynomials. 

Jointly  with  S.  Gottlieb,  we  have  further  explored  a  class  of  high  order  TVD 
(total  variation  diminishing)  Runge-Kutta  time  discretization  initialized  by  Shu  and 
Osher,  suitable  for  solving  hyperbolic  conservation  laws  with  stable  spatial  discretiza¬ 
tions  [18].  We  illustrate  with  numerical  examples  that  non- TVD  but  linearly  stable 
Runge-Kutta  time  discretization  can  generate  oscillations  even  for  TVD  (total  vari- 
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ation  diminishing)  spatial  discretization,  verifying  the  claim  that  TVD  Runge-Kutta 
methods  are  important  for  such  applications.  We  then  explore  the  issue  of  optimal 
TVD  Runge-Kutta  methods  for  second,  third  and  fourth  order,  and  for  low  storage 
Runge-Kutta  methods. 

We  have  been  continuing  our  investigation  on  adapting  shock  capturing  algorithms 
for  semiconductor  device  simulations.  Jointly  with  Jerome  [3],  we  have  studied  the 
effect  of  the  common  practice  of  neglecting  the  convective  terms  (inertial  approxi¬ 
mation)  in  the  hydrodynamic  model  in  device  simulations.  We  find  out  that  this 
inertial  approximation,  though  very  convenient  for  the  application  of  exponential  fit¬ 
ting  type  methods,  is  nevertheless  not  valid  near  the  diode  junctions,  resulting  in 
significant  error  in  the  result.  We  have  also  studied  the  characteristic  modes  of  some 
modified  hydrodynamic  models  and  energy  transport  models,  the  purpose  being  to 
make  sure  that  such  modifications  do  not  change  the  mathematical  properties  of  the 
models  dramatically.  Jointly  with  Chen,  Cockburn  and  Jerome,  we  have  extended 
the  discontinuous  Galerkin  method  to  the  hydrodynamic  model  of  semiconductor 
device  simulations  [6].  Jointly  with  Jerome,  we  have  studied  the  response  of  the 
hydrodynamic  model  to  heat  conduction,  mobility,  and  relaxation  expressions  [7]. 
Also  jointly  with  Jerome,  we  have  performed  analysis  and  simulation  for  the  energy 
transport  models  for  semiconductors  [9]. 

Jointly  with  G.  Chen,  J.  Jerome  and  D.  Wong  of  Northwestern  University,  we 
have  introduced  a  novel  two  carrier  (electro)  hydrodynamic  model  for  semiconductor 
device  simulations,  which  incorporates  higher  dimensional  geometric  effects  into  a  one 
dimensional  model[19],  [20].  A  rigorous  mathematical  analysis  is  carried  out  for  the 
evolution  system  in  the  case  of  piezotropic  flow,  including  realistic  carrier  coupling. 
The  proofs  are  constructive  in  nature,  making  use  of  generalized  Godunov  schemes 
with  a  novel  fractional  step,  steady-state  component,  and  compensated  compactness. 
Two  important  applications  are  studied.  We  simulate:  (1)  the  GaAs  device  in  the 
notched  oscillator  circuit;  and,  (2)  a  MESFET  channel,  and  its  steady-state  symme¬ 
tries.  The  first  of  these  applications  is  the  well  known  Gunn  oscillator,  and  we  are 
able  to  replicate  Monte-Carlo  simulations,  based  upon  the  Boltzmann  equation.  For 
the  second  application,  we  observe  the  effect  of  a  symmetry  breaking  parameter,  the 
potential  bias  on  the  drain.  In  high  order  essentially  non-oscillatory  (ENO)  finite 
difference  schemes, 

Jointly  with  Cercignani,  Gamba  and  Jerome,  we  have  been  investigating  a  new 
high  field  model  for  the  semiconductor  devices  [23],  [24].  It  is  our  hope  that  this  new 
model,  coupled  with  a  domain  decomposition  technique,  will  give  more  accurate  de¬ 
scription  of  the  small  devices  currently  being  used  with  relatively  small  computational 
effort. 

Jointly  with  Chen,  Eisenberg  and  Jerome,  we  have  explored  a  hydrodynamic 
model  of  temperature  change  in  open  ionic  channels  [8].  This  is  an  important  problem 
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in  computational  biology. 
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